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ABSTRACT 

We present multiple ultra-high resolution cosmological hydrodynamic simulations of M* ~ 
104“®-3]\/[^ dwarf galaxies that form within two Mvir = 10®-3“40]y[^ matter halo ini¬ 
tial conditions. Our simulations rely on the FIRE implementation of star formation feedback 
and were run with high enough force and mass resolution to directly resolve structure on the 
^ 200 pc scales. The resultant galaxies sit on the M* vs. Mvir relation required to match 
the Local Group stellar mass function via abundance matching. They have bursty star forma¬ 
tion histories and also form with half-light radii and metallicities that broadly match those 
observed for local dwarfs at the same stellar mass. We demonstrate that it is possible to cre¬ 
ate a large 1 kpc) constant-density dark matter core in a cosmological simulation of an 
M* ~ 1O®’3]V[0 dwarf galaxy within a typical Mvir = IO^^Mq halo - precisely the scale of 
interest for resolving the Too Big to Fail problem. However, these large cores are not ubiq¬ 
uitous and appear to correlate closely with the star formation histories of the dwarfs: dark 
matter cores are largest in systems that form their stars late (z < 2), after the early epoch 
of cusp building mergers has ended. Our M,^ ~ IO^Mq dwarf retains a cuspy dark matter 
halo density profile that matches that of a dark-matter only run of the same system. Though 
ancient, most of the stars in our ultra-faint form after reionization; the UV field acts mainly to 
suppress fresh gas accretion, not to boil away gas that is already present in the proto-dwarf. 

Key words: galaxies: formation — galaxies: evolution — galaxies: dwarf — cosmology: 
theory — methods: numerical 


1 INTRODUCTION 

Many of the most pressing problems associated with the standard 
LCDM paradigm concern the faintest M* ~ IO^Mq dwarf galax¬ 
ies and the dark matter halos that have the right abundance to host 


them: Mvir — 10 “Mq (Garrison-Kimmel et al. 


2014 


Brook et al. 


|2014| >. If LCDM is correct, then the dark matter halos hosting these 
dwarfs must be extremely inefficient at converting baryons into 
stars jKlypin et al.|199^ and they also must be significantly less 
dense in their centers than predicted in dissipationless LCDM sim- 


ulations (Boylan-Kolchin et al. 

2011 

2012 

Ferrero et al.||2012[ 

Garrison-Kimmel et al.||2014 

blleruc 

et a 

20T^ Klypin et al'| 

2014| Papastergis et al.|2015 1 . This latter issue (known as the Too 


E-mail: onorbe@mpia, de 


Big to Fail problem) may be related to indications that dwarf galax¬ 
ies reside within dark matter halos that have cored density profiles 
rather than the cuspy NFW-like profiles predicted in CDM simula- 


tions (|Flores & Primack|1994 

Kuzio de Naray et al.|20081 de Blok 

et al.|2008l|0h et al.|2008l|Wa 

cer & Penarrubia|201 lllSalucci et al. 

2012||Amorisco et al.|2014||Ogiya & Burkert|20151 but 
et al. 2014). 

While some authors have taken these discrepanc 
tivation to explore non-standard dark matter models i 

»ee Strigari 

les as mo- 

Maccio & 

Fontanot||20101 Vogelsberger et al.| 2012| |Rocha et al.| 

20131 |Ho- 


that it may be possible to naturally resolve them through a better 
understanding of star formation and feedback in low-mass galaxies. 
Specifically, the inefficiency of dwarf galaxy formation is believed 
to be driven by supemovae feedback and the effects of an ioniz- 
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2 J. Onorbe et al. 


ing background jPekel & Silk|[T98^ |Bullock et al.||20Q0) . Like¬ 
wise, dark matter halos may be transformed from cusps into cores 
if enough energy can be injected into the orbits of dark matter par¬ 
ticles during rapid starburst events ([Nayarro^t alJ]996 J 


|et al.|201^|Pontze n & Govemato|2012 Ogiya & Mori|2014 1 . As 

pointed out by |Penarrubia et al.| ( |2012[ l, these two requirements are 
at odds with each other: the need to lower the efficiency of star for¬ 
mation means that there will be less supernovae energy available to 
create dark matter cores. Solving the two problems simultaneously 
therefore represents a significant theoretical challenge jGarrison^ 
[Kimmel et al.|20T^ . 

Reproducing even the broad-brush properties of dwarfs in a 
cosmological framework, regardless of their internal structure, has 
heen historically challenging. At these scales, the relationship be¬ 
tween stellar mass and halo mass derived from local galaxy counts 
( [Garrison-Kimmel et al.|2014l|Brook et al.|2014^ implies a suppres¬ 
sion of galaxy formation hy a factor of about 1000. While it is gen¬ 
erally believed that stellar feedback is the main agent responsible 
for this suppression, actually getting a physically realistic model 
of the relevant processes to manifest these expectations has proven 
difficult. 

The past several years have proven fruitful in this regard, with 
many published studies achieving substantial suppression in the 


conversion of baryons to stars on the scale of dwarf galaxy halos 

([Governato et al.[[2010 Sawala et al. 

2011[[Simpson et al.|20131 

Munshi et al.|2013| Governato et al. 

2015| [Trujillo-Gomez et al. 

2015[l. As we show below, however. 

many of these studies have 


not quite reached the level of suppression that seems to he required 
hy local galaxy counts. Moreover, whether or not these feedback 
models also match the different observed scaling relations for these 
systems ( |Wolf et al.|20T0l [Kirby et al.|2013[[Cdlins et al.|2014| ) is 
still not clear. Reproducing both the correct stellar mass and struc¬ 
tural properties has proven to be an even more difficult challenge 
([Sales et al.|2010^. The observed stellar metallicity - stellar mass 


tight correlation (Gallazzi et al.|2005| [Kirhy et al.|2013) can also 


put very important constraints on the feedback models and how 
these are implemented. 

As for the question of feedback-driven core formation, much 
remains debated. Some of the most successful simulations at pro¬ 
ducing cores in dwarf galaxies have suggested a transition ma ss be- 
low M* ~ 10^ Mq where core formation becomes difficult jGov- 


emato et al. 


Cintio et al. 


201^L Using a slightly different set of simulations, jPi 


1 2014|( find similar results, and suggest that the cusp- 


core transition should be most effective when the ratio of stellar 
mass to dark matter halo mass relatively high, in massive dwarfs 
with M* ~ 10® M© and Mvir ~ 3 x Mq. Importantly, 

they also find that cuspy profiles are retained for the M* ~ IO^Mq 
dwarfs of concern (residing in Mvir = IO^^Mq halos) although 
resolution may have been an issue in these cases. At some mass 
scale, galaxy formation may become effectively stochastic (e.g., 
Boylan-Kolchin et al.|2011^. Recent work by jSawala et al.|j2bl5| 


2014 1 , however, suggests that the scale at which stochasticity be¬ 
comes important is somewhat lower (Mvir ~ 10® - 1O® ®M0). 

Though the results of jDi Cintio et al.H2014^ and |Governato| 
jet al.j \2Q\2) agree reasonably well, a different set of high reso¬ 
lution simulations with a simpler implementation of stellar feed¬ 
back have not produced cores in dwarf galaxy halos at any mass 
( jVogelsberger et al.|2014| l, even though a number of other observ¬ 
ables are well matched. The absence of cores produced by stellar 
feedback in these simulations could be due to the fact that their 
sub-grid ISM and star formation model leads to star formation his¬ 
tories that are (likely) artificially smoothed in time, compared to 


the bursty star formation histories found in more explicit models 
([Hopkins et al.||2014[|Muratov et al.||2015^ . Conversely, [Trujill^ 
[Gomez et al.[j2015[ l found that radiation pressure from massive 
stars was the most important source of core formation in their sim¬ 
ulations, not thermal feedback from supernova, which has been the 
primary mode used by other groups that have produced cores. More 
generally, models for feedback that have been used up until now 
have been sub-grid and necessitated ad-hoc approximations, such 
as turning off cooling for material heated by SNe. As such it is 
not clear whether the feedback we actually expect from stellar evo¬ 
lution models is capable of producing large cores, or whether the 
mass-limit for core formation is robust. 

In this paper, we attempt to minimize the freedom of sub¬ 
grid galaxy formation models and to incorporate as many impor¬ 
tant physical processes in a manner that is as realistic as possible at 
present in order to understand if and how star formation affects the 
gravitational potential wells of dwarf dark matter halos. To these 
ends, we have conducted a series of high resolution cosmologi¬ 
cal hydrodynamical simulations of two dwarf halos using the code 
presented in [Hopkins et al.[ ( [2014[ (. In this work, we showed that 
this implementation of stellar feedback successfully reproduces the 
observationally-inferred relationship between the stellar mass-dark 
matter halo mass (M*-Mhaio) and star formation histories of galax¬ 
ies at all redshifts where observational constraints are currently 
available. [Faucher-Giguere et al.[ ( [2015t recently showed that it also 
replicates the neuhal hydrogen content of high-redshift halos. 

To our knowledge, the set of simulation presented here include 
the current highest resolution simulation of this type with an ex¬ 
plicit implementation of feedback yet achieved. This not only fa¬ 
cilitates a more accurate treatment of astrophysical processes but 
is also crucial in the context of dwarfs as dark matter probes. The 
dwarfs of concern have half-light radii of ~ 500 pc, and thus any 
dark matter core of relevance needs to be dynamically resolved at 
this scale. According to well-documented convergence test studies 
( [Power et al.|2003] l, many previous simulations that have reported 
core formation on this scale were quite poorly resolved, some at 
only ~ 2 — 3 softening lengths. In what follows we make every 
effort to clarify our resolution limitations. 

The paper is organized as follows. In Section we describe 
the computational methods which we have used and our choice of 
initial conditions. We present the results of our simulations in Sec¬ 
tion]^ We pay closer attention to the matter content of our simu¬ 
lated dwarfs, and the possible formation of cores, in Section|^ We 
conclude with a summary where we discuss the achievements and 
shortcomings of the simulations in Section]^ 


2 SIMULATIONS 

We have run a series of multimass cosmological hydrodynamical 
simulations ( [Porter[19M|[Katz & White[1993^ following the forma¬ 
tion and evolution of structure in the ACDM model of two dwarf 
galaxy halos. Each simulation is a cosmological zoom-in that in¬ 
cludes high-resolution gas and dark matter for the flow converging 
region that generates the main object. The rest of the simulation box 
is sampled by low-resolution dark matter particles that account for 
tidal forces. The cosmological model adopted throughout this paper 
is based on cosmic microwave background results ([Komatsu et al.j 
|20TT] i: erg = 0.801, Ha = 0.734, = 0.266, Hb = 0.0449, 

Us = 0.963 and h — 0.71. 

To generate the cosmological initial conditions we made use 
of MUSIC, an OPENMP parallel algorithm to generate multi-scale 


© 0000 RAS, MNRAS 000, 000-000 














































































Forged in FIRE 3 



Figure 1. From left to right, visualizations of the gas density, gas temperature and gas metallicity for the DwarLearly run at 2 : = 2.3. All panels show the 
same thin slice along the 2 -axis centered at the main halo. The signatures of a recent stellar burst episode are clear in all of them. 


initial conditions with multiple levels of refinements for cosmo¬ 
logical “zoom” simulations (MUSIC Hahn^&^bel j^lT} and we 
followed the method outlined in |Onorbe et al.| ( |2014| ). To select 
our dwarf candidates we first run a medium-resolution dark-matter 
only cosmological simulation using GADGET-2 jSpringef||2005^ 
with a cubic volume of 7 Mpc on a side with particle mass nip = 
9.7 X and Plummer equivalent force softening length of 

176 pc. To be able to study the main statistical properties of dwarf 
galaxy halos we also run a bigger dark-matter only simulation of 35 
Mpc on a side with particle mass trip = 1.2 x IO^Mq and Plummer 
equivalent force softening length of 563 pc. In this work we present 
simulations of two dwarf galaxy halos, one with a virial mass of 
Mvir = 3.2 X 10®Mq and the other with Mvir = 9.2 x 10®M(r[^ 
Based on our analysis of the 35 Mpc simulation, we have chosen 
our dwarf candidates to lie as close as possible to the mean values 
of spin, concentration and halo formation time for its mass while 
still having a small Lagrangian volume (see |Onorbe et al.|2014^ . 
The specific values of these parameters for our two halos can be 
found in TableWe point to Appendix]^ for a more detailed de¬ 
scription of these parameters and how they compare with a sample 
of halos in the same mass bin. 

To check the convergence of our results we have run two 
resolution levels for our simulations: in our low-resolution hy- 
drodynamical testing runs we use a dark matter particle mass of 
1.01 X 10"^ Mq and a particle gas mass of 2.04 x 10^ M© (the 
mass resolution for the collisionless run is therefore 1.22 x lO'* 
M©). The high resolution runs used a dark matter particle mass of 
1.26 X 10® M© and a gas particle mass of 254 M© (the particle res¬ 
olution for the collisionless run is therefore 1.5 x 10® M©). None 
of the high resolution regions of the simulations presented in this 
work are contaminated by low resolution particles at any redshift 

within 1.6 virial radii. _ 

The simulations presented in this paper use GIZMC[^ (Hop 


|ldEi||20l4l l, run in P-SPH mode which include physical models 
for star formation and stellar feedback presented in|Hopkins et al.| 
( |20T4) . Two of the runs presented here Ultrafaint and DwarUearly 
were also presented in [Hopkins et al.|f2014| l (m09 and mlO respec¬ 
tively). We summarize their properties below, but readers interested 
in further details (including resolution studies and a range of tests 


of the specific numerical methodology) should see|Hopkins et al.| 

( |M2l[2m3l|M4t . 

For the halo identification in the simulation we have used the 
public code Amiga Halo Finder (AHF |Knollmann & Knebe|2009^ , 
an MPI parallel code for finding gravitationally bound structures in 
simulations of cosmic structure. Results presented in this work use 
a highest density peak-l-sigma-clipping method to find the center. 
We have also tested different centering algorithms to confirm that 
our results do not depend on which method was usecj^ 


2.1 Numerical Methods 

The P-SPH method adopts the Lagrangian “pressure-entropy” for¬ 
mulation of the SPH equations developed in jHopkins|2013^ ; this 
eliminates the major differences between SPH, moving mesh, and 
grid (adaptive mesh) codes, and resolves the well-known issues 
with fluid mixing instabilities in previously-used forms of SPH 
(e.g. |Agertz et al.|2007[|Sijacki et al.|2012^ . P-SPH also manifestly 
conserves momentum, energy, angular momentum, and entropy. 
The gravity solver is a heavily modified version of the GADGET- 
3 ( |Springel|2005^ hybrid tree-particle mesh (Tree-PM) method; but 
GIZMO also includes substantial improvements in the artificial vis¬ 
cosity, entropy diffusion, adaptive timestepping, smoothing kernel, 
and gravitational softening algorithm, as compared to the “previous 
generation” of SPH codes. These are all described in detail in |Hop-| 
|kins et aL| ( |2014) ; |Hopkins| j2014^ . In particular, in “traditional” 
GADGET, softenings are not adaptive, and pairwise interactions 
are simply smoothed by the larger of the two particle softenings. 
We have also modified the softening kernel as described therein to 
represent the exact solution for the potential of the SPH smoothing 
kernel. Therefore our “standard” simulations use adaptive gravita¬ 
tional softening lengths for gas which minimum is a factor ~ 10 
smaller than the fixed dark matter gravitational softening lengths. 
In order to test this approach we have also run the same initial con¬ 
ditions using identical softenings for both the baryonic and dark 
matter particles (close to the higher dark matter default value). We 
labeled these runs according to the late star formation history of the 
high resolution runs (see Table[T]and the discussion below for more 
details). 

In our simulations, gas follows an ionized-tatomic-l-molecular 


^ Unless otherwise stated, in this paper we define the virial overdensity 
using the spherical top hat collapse approximation by |Bryan & Norman| 

://www.tapir.caltech.edu/ phopkins/Site/GIZMO 



® A simple center-of-mass algorithm was the only method that we found 
not able to track the center of our systems with the accuracy required for 
this work. 
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cooling curve from 10 — 10^° K, including metallicity-dependent 
fine-structure and molecular cooling at low temperatures, and high- 
temperature (> 10^ K) metal-line cooling followed species-hy- 
species for 11 separately tracked species. At all times, the ap¬ 
propriate ionization states and cooling rates are tabulated from 
a compilation of CLOUDY runs, including the effect of a uni- 
for m but redshift-dependent photo-io nizing background computed 
in I Faucher-Giguere et al. 2009|'^ together with local sources 
of photo-ionizing and photo-electric heating. Self-shielding is ac¬ 
counted for with a local Sobolev/Jeans-length approximation (in¬ 
tegrating the local density at a given particle out to a Jeans length 
to determine a surface density E, then attenuating the background 
seen at that point by exp(«:yE)). 

Star formation is allowed only in dense, molecular, self- 
gravitating regions above n > ricrit (ricrit = 100 cm“® for 
our high-resolution simulation^. This threshold is much higher 
than that adopted in most “zoom-in” simulations of galaxy forma¬ 
tion (the high value allows us to capture highly clustered star for¬ 
mation). We follow IKrumholz & Gnedin| j201 1^ to calculate the 
molecular fraction /h 2 in dense gas as a function of local column 
density and metallicity, and allow SF only from molecular gas. We 
also follow [Hopkins et al.H2013} and restrict star formation to gas 
which is locally self-gravitating, i.e. has a = 5v^Sr/G'm^^s{< 
Sr) < 1 on the smallest available scale (Sr being our force soften¬ 
ing or smoothing length). This forms stars at a rate p* = pmoi/ts 
(i.e. 100% efficiency per free-fall time); so that the galaxy and even 
kpc-scale star formation efficiency is not set by hand, but regulated 
by feedback (typically at much lower values). 

Feedback from stellar evolution is modeled by implementing 
energy, momentum, mass, and metal return from radiation, super¬ 
novae, stellar winds, and photoionization. Every star particle is 
treated as a single stellar population, with a known age, metal¬ 
licity, and mass. Then all feedback quantities (the stellar luminos¬ 
ity, spectral shape, SNe rates, stellar wind mechanical luminosities, 
metal yields, etc.) are tabulated as a function of time directly from 
STAREURST99 stellar population synthesis model jLeitherer et al.| 
|1999| l, assuming a |Kroupa]j2002^ IMF. Details on the implemen¬ 
tation of each of these physical processes in our simulations can 
be found in ( [Hopkins et al.|20l4) . No black hole physics has been 
considered in these simulations. 

Despite taking all our inputs directly from stellar population 
models, there are some ambiguities in how we implement them. 
For example, when we deposit mass, momentum, and energy to 
particles within the SPH kernel, we can do so according to a mass¬ 
weighting or volume-weighting scheme. We have experimented 
with both, and we refer to these options as Feed-M and Feed-V, 
respectively. 

We stress that the systematic differences due to these (and 
other similar) purely numerical choices (see Appendix A ofj Hop-j 
[kins et al.|2014| l are relatively small for integrated quantities like the 
stellar mass. However, since the dynamics of galaxies and star for¬ 
mation are chaotic, a small perturbation can make a non-negligible 
difference to the shape of the star formation history. These essen¬ 
tially stochastic variations will provide a useful means for us to ex¬ 
amine the role of different star formation histories in shaping cores. 

We have found that the main global parameters describing the 
dwarf galaxies are quite robust regardless of resolution, softening 


^ Publicly available at http://galaxies.northwestem.edu/uvb 
® The simulations get to a maximum density of ~ lO* cm“® 


and other minor changes in the code. See Appendix [^ for a full 
discussion on the convergence of our results. 


2.2 Sample Summary 

A summary of all the relevant parameters used in the ultra high res¬ 
olution runs presented in this work is shown in Table [T] along with 
the naming conventions we have adopted. In this work we present a 
total of six high resolution simulations of two dwarf galaxy halos, 
one with a virial mass of Mvir = 3.2 x IO^Mq and the other with 
Mvir = 9.21 X 10® Mq (as measured in the high-resolution col¬ 
lisionless simulations). For the more massive halo we present here 
four runs, a high resolution collisionless run (DwarLdm) and a total 
of three different hydrodynamical runs which include two feedback 
implementation tests and the softening test mentioned above. 

We have named the three hydrodynamical Dwarf simulations 
based on their star formation histories (see section 3.2 below). 
The run we call “Dwarf_early” shows most of its star formation 
at early times and corresponds with the feedback method Feed- 
V. The run we call “DwarfJate” uses feedback method Feed-M 
and shows a more significant star formation rate at low redshifts. 
The “Dwarfjniddle” run is the softening test which uses feedback 
method “Feed-M” and its star formation rate history stands just be¬ 
tween the two. Simulations of the same dwarf using the “Meshless 
Finite Mass” method implemented in GIZMO ( |Hopkins|2014[ ( and 
the feedback Feed-V) method produce results very similar to the 
“Dwarf_early” run presented here (Fitts et al., in preparation). 

For the smaller halo we have run the same number of simula¬ 
tions as we have for the larger one, but their results were so similar 
that we present only one hydrodynamic run (Ultrafaint, which uses 
Feed-V) and one collisionless run (Ultrafaint_dm). The hydrody¬ 
namical runs DwarLearly and Ultrafaint were already presented in 
the first FIRE paper ( [Hopkins et al.|2014[ . 

As pointed out above, we have also checked the convergence 
of these results with resolution by running all these setups also at 
a lower resolution level. We discuss these runs in detail in Ap¬ 
pendix [^ We have also run many more (~ 50) simulations at 
this lower resolution level of these halos to test other purely nu¬ 
merical issues and the effects of adding/removing each feedback 
mechanism in turn. Some of these are summarized in [Hopkins et al.[ 
( [2014[ . We will not discuss them further in this paper because they 
are either not instructive for the study of this work because the in¬ 
cluded physics is not complete or because there is no change in the 
results. Even given excellent force and mass resolution, the time- 
step criterion used in simulations is always a concern if many, many 
orbits of N-body particles must be followed (as in the halo cen¬ 
ters of the systems studied here). These can artificially deteriorate a 
central cusp, if an insufficiently stringent timestep criterion and/or 
error tolerance for the long-range force computations is used. We 
have therefore re-run a subset of our low-resolution runs, making 
the timestep criterion a factor of ~ 30, and force error tolerance 
a factor of ~ 100 times more strict than our default choices. This 
amounts to taking < 100 year timesteps, with a tree force accuracy 
a factor ~ 1000 stricter than used in [Governato et al.] ( [2012[ (, and a 
factor ~ 100 stricter than was found to give good convergence in 
idealized comparisons of dark matter zoom-in simulations in ( [Kim[ 
[et al.[2014| l. Given our very strict default tolerances, this gave well- 
converged results. 

Figure shows visualizations of the gas density (left panel), 
gas temperature (middle panel) and gas metallicity (right panel) for 
the DwarLearly run at 2 = 2.3. All panels show the same thin 
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Figure 2. Galaxy stellar mass-halo mass relation at z = 0 for the simula¬ 
tions presented in this work. Filled red triangle stands for the Ultrafaint run. 
Empty red square, filled red pentagon and filled red diamond stand for the 
Dwaif.early, DwarLmiddle and Dwarfjate runs respectively. Circles stand 
for other published simulations in the literature using sub-grid stellar feed¬ 
back models. Empty symbols stand for runs in which a dark matter core was 
found according to the typical definition used in previous works: a slope 
a > —0.9 for the dark matter density profile between 1% - 2% of the virial 
radius. The two black lines show abundance matching relations derived by 
Brook et al. and Garrison-Kimmel et al. using galaxy counts within the Lo¬ 
cal Volume. These relations should be complete to M* ~ 10 ®Mq (much 
deeper than other published abundance matching relations and therefore the 
most relevant for this comparison). To make all simulation data and abun¬ 
dance matching results consistent between each other, we have plotted the 
total virial masses of the dark matter only runs defined as Avir = 200pcrit ■ 
A small correction has been applied whenever these values were not avail¬ 
able in the literature (with 10% effects). In order to facilitate direct com¬ 
parison with other results in the literature, the stellar mass plotted is all of 
the stellar mass inside the virial radius. Using a smaller radius makes only 
a small difference for our simulated dwarf galaxies. 


slice along the z-axis centered at the main halo. The signatures of 
a recent SN episode are clear in all of them. 


mass range presented here, this is particularly impressive, as the in¬ 
tegrated stellar mass is suppressed by factors of ~ 1000 relative to 
the Universal baryon fraction (upper solid gray line). 

The smaller points in Figurej^show results from previous hy- 
drodynamical simulations of dwarf galaxies (Governato et al.|! 


Sawala et al.|2011[|Simpson et aL|2013[ |Munshi et al.|2013( ^en| 


et al.|2014| Trujillo-Gomez et al.|2015 i. The open points are those 

that have reported at least mild flattening of the central dark matter 
cusp in response to feedback effects. We note that all of those open 
points are associated with systems that have formed a fair number 
of stars, with M* > 7 x 10® Mq - more massive than the systems 
of concern for the Too Big to Fail Problem. As we discuss below, 
one of our runs (Dwarfdate, open square) produces a large core 
while forming significantly fewer stars. 

The upper left panel of Figure]^ shows galaxy size, measured 
as the half-stellar mass radius, versus the total stellar mass of the 
galaxy for our simulated galaxies (red points). The observed stellar 
size-mass relation seen for Milky Way satellites (green) and Local 
Group field dwarfs (yellow) are shown as data points (taken from 
|Wolf et al.|2010HKirby et al.|2013[|2014[ who also compile data 
from the literature). The upper right panel of Figurej^shows the to¬ 
tal mass within the half-stellar mass radius vs. the half-stellar mass 
radius, again for our simulated galaxies compared to local galax- 
ieQ Finally, the bottom panel shows the stellar-mass metallicity 
relation. Overall, the simulated galaxies are in good agreement with 
sizes, metallicities, and total masses seen for galaxies of their stellar 
mass in the local universe. For example, the Dwarf date and middle 
runs show a good agreement with Fornax. Dwarf_early’s size and 
mass are close to Ursa Minor. 

Table [T] shows that DwarLearly is much more dark-matter- 
dominated within 500 pc than Dwarfdate. This also holds when 
we look at the half-stellar mass radius of each galaxy, instead 
of at a fixed physical value. Within this radius, DwarLearly has 
Mtot/M* « 44 and Mtoti/Mbaryons « 11. By mass, stars are 
subdominant to gas within the half-light radius by a factor of 2.8 
for this dwarf. Within the half-mass stellar radius, Dwarfdate has 
Mtot/M* Ri 87 and Mtot/Mbaryons « 5.4 owing to a large reser¬ 
voir of gas within the stellar half-mass radius (Mgaa = 16M* 
within r*/ 2 )- 


3 RESULTS 

3.1 Basic Properties at z = 0 

Table [T] presents some relevant parameters describing the proper¬ 
ties of each simulation presented in this work that will allow an 
immediate comparison with previous simulations and observations 
of dwarf galaxies. 

Of particular interest is the resultant stellar mass in each 
dwarf. Figure|^presents the stellar mass - halo mass relation for the 
four hydrodynamical runs described above (large red points) com¬ 
pared to the most recent estimates for this relation from abundance¬ 
matching exercises in the Local Group jGarrison-Kimmel et al.| 
|2014[[Brook et al.|2014| black solid and dashed lines respectively). 
The known sources of stellar feedback we include, with no adjust¬ 
ment, automatically produces galaxy stellar masses that are consis¬ 
tent with those required to match local galaxy count^ For the halo 

® Abundance matching results below ~ IO^Mq stellar mass are extrapo- 


3.2 Star Formation Histories 

While the simulated Dwarfs (early, middle, and late) all show sim¬ 
ilar z = 0 stellar masses, they arrived at those final states via dif¬ 
ferent paths. The Ultrafaint run, on the other hand, ends up with a 
stellar mass some two orders of magnitude smaller than any of the 
Dwarf runs, though it resides within a halo that is only ~ 3 times 
less massive. In this subsection, we explore these differences by 
examining the star formation histories in some detail. 

In Figure]^ we present the star formation rates (left panel) and 
the normalized cumulative star formation histories (right panel) of 
all four high resolution hydro runs. The Ultrafaint simulation (or¬ 
ange line) forms all of its stars before z ~ 2.5. The galaxy shuts 
down at this redshift because of two main effects; 1) the UV back¬ 
ground prevents fresh gas accretion after z ~ 6 and 2) stellar feed¬ 
back acts to self-quench the system after the ionizing background 


lation as observations throughout the Local Group are not complete below 
this limit 

^ When we do this comparison with observations we are assuming that the 
half-stellar mass radius is equivalent to the half-light radius. 
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Figure 3. All panels show the hydrodynamical runs presented in this work compared with different observations aX z = 0. Empty red square, filled red 
pentagon and filled red diamond stand for the DwarLearly, DwarLmiddle and DwarfJate runs respectively. Upper left panel: Stellar mass versus effective 
stellar mass radius. Circles stand for observational data from jWolf et al.|2010||Kirby et al.|2013||2014) . Upper right panel: The galaxy size (effective stellar 
mass radius) versus mass (total mass inside this radius) relation for the simulations presented in this work. In order to compare it with observations, we have 
plotted the data derived by |Wolf et al.HlOTO) (green circles) and |Kirby et al.|j2014j ^(orange circles). Bottom panel: Stellar mass - stellar metallicity relation 
for the simulations presented in this work. Circles stand for the data compiled by|Kirby et al.|j2013^ for nearby dwarfs. See text for more details. 


turns on. The remaining cold gas found at 2 ; = 0 is not able to 
reach high enough densities to generate stars. Also notice that even 
though the UV background starts acting at high redshift on the gas 
particles, there is still active star formation for about one billion 
years after this redshift. 

This is not unexpected: previous lower-resolution simula¬ 
tions predicted that reionization-induced UV heating is not strong 
enough to remove all of the gas from dwarf-sized halos jHoeft et al.| 
|2006| l. However, these simulations were not able to resolve the star 
formation histories of these galaxies, so it was not clear if the re¬ 
maining gas would be able to form stars. As the UV background 
effectiveness depends on the density of the gas, cold and dense gas 
is not affected. Moreover, this more efficient star formation period 
seems crucial in order to match the stellar metallicity ratios ob¬ 
served for low-mass dwarf galaxies ( |Kirby et al.|2013^ . 

Figure|^also shows the star formation rates of the three Dwarf 
runs. Dwarf_early forms more than half of its stars prior to 2 = 2.5 
while DwarfJate maintains a fairly substantial star formation rate 
down to 2 = 0. The Dwarf_middle star formation rate history 


stands just between the two. All of the runs show bursty star forma¬ 
tion histories on ~ 100 Myr timescales. 

In all three of the Dwarf runs, the star formation histories 
show two different phases. At the highest redshifts (2 > 3), total 
dark halo and stellar masses both grow efficiently (albeit with some 
offset). This is the “rapid assembly” phase ( [Wechsler et al.|2002| >, 
before/during reionization, in which feedback, while able to eject 
some gas from the galaxy and provide some overall suppression 
and variability of the star formation, does not appear to dominate 
the gas dynamics (the central potential and mass of the halo grow 
on timescales comparable to the galaxy dynamical time). But from 
2 ~ 3 onward, halo accretion rates slow down and feedback acts 
strongly. From this point on, there appears to be a steady-state SFR 
that can be considered constant with time when averaged over a 
Gyr scale (a bursty behavior emerges when smaller time bins are 
used). In this phase, the galaxy is able to cycle new material into 
a fountain and so maintain equilibrium. This “quasi-equilibrium” 
SFR scales with the central potential of the galaxy (see |Hopkins| 
|et al.|2012| ), as traced by quantities such as the central halo den¬ 
sity or Vmax (the maximum circular velocity), not the halo mass 
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Parameter 

Ultrafaint_dm 

Ultrafaint^ 

Dwarf _dm 

Dwarf_late 

Dwarf .middle 

Dwarf.early^ 


(Collisionless) 

(Hydro: Feed-V) 

(Collisionless) 

(Hydro: Feed-M) 

(Hydro: Feed-M-soft) 

(Hydro: Feed-V) 

(Mq) 

1.5 X 10® 

1.26 X 10® 

1.5 X 10® 

1.26 X 10® 

1.26 X 10® 

1.26 X 10® 

2) Edm (pc) 

28 

28 

35 

35 

25 

28 

3) (Mq) 

- 

2.54 X 10® 

- 

2.54 X 10® 

2.54 X 10® 

2.54 X 10® 

4) (pc) 

- 

1.0 

- 

2.0 

25 

2.8 

5) Mvir (Mq) 

3.2 X 10® 

2.5 X 10® 

9.2 X 10® 

7.6 X 10® 

7.7 X 10® 

7.7 X 10® 

6) Vmax (km/s) 

26 

22 

37 

33 

33 

32.5 

7) Tyir (kpc) 

38 

35 

54 

51 

51 

51 

8) A 

0.031 

- 

0.0350 

- 

- 

- 

9) Vmax!^vir 

1.35 

- 

1.38 

- 

- 

- 

10) *50 (Gyr) 

1.43 

- 

1.84 

- 

- 

- 

il) /bar X (S2ni/flb) 

- 

0.024 

- 

0.093 

0.074 

0.056 

12) /gas 

- 

0.0049 

- 

0.018 

0.014 

0.011 

13)/* 

- 

0.00002 

- 

0.0004 

0.0003 

0.0003 

14) M* (Mq) 

- 

2.1 X 10"^ 

- 

2.8 X 10® 

2.7 X 10® 

2.2 X 10® 

I5)r*^2 (P^) 

- 

340 

- 

1100 

830 

550 

16) [Fe/H] 

- 

-2.8 

- 

-2.0 

-2.0 

-1.9 

17) M^^t, (Mq) 

2.2 X 10 '^ 

1.7 X 10® 

3.4 X 10® 

0.75 X 10® 

1.3 X 10® 

1.9 X 10® 

18) Mfg^k (Mq) 

1.6 X 10 '^ 

1.7 X 10® 

2.6 X 10® 

0.43 X 10® 

0.80 X 10® 

1.6 X 10® 

19) Mb- (Mq) 

- 

1.2 X 10'^ 

- 

3.2 X 10® 

4.7 X 10® 

2.2 X 10® 

20) M|“ (Mq) 

- 

6.9 X 10® 

- 

3.1 X 10® 

4.5 X 10® 

1.6 X 10® 

21) M*oo (Mq) 

- 

1.9 X 10'^ 

- 

5.4 X 10"^ 

2.4 X 10® 

6.3 X 10® 


Table 1. Simulations data. First column stand for the different parameters studied for each simulation. In Columns 2-7 results at 2 = 0 for the simulations 
presented in this work are shown. Row 1: dark matter particle mass in the high resolution region in solar masses. Row 2: fixed gravitational softening used 
for the dark matter particles in physical pai'secs. Row 3: baryon particle mass in the high resolution region in solar masses. Row 4: minimum baryonic force 
softening in parsecs (minimum SPH smoothing lengths are comparable or smaller). Recall that force softenings are adaptive (mass resolution is fixed). Row 
5: virial mass in solar masses defined at the overdensity at which the spherical top hat model predicts virialization jBryan & Norman|1998^ . Row 6: maximum 
circular velocity in km/s. Row 7: virial radius in kiloparsecs. Row 8: halo spin (Bullock et al.|20QT} definition. See Appendix|A|for more details. Row 9: halo 
concentration. See Appendix[^for more details. Row 10: halo formation time. See Appendix|A|fm' more details. Row 11: virial baryon fraction, i.e., baryon 
mass inside the virial radius over the virial mass, divided by the cosmic baryon fraction. Row 12: virial gas fraction, i.e., gas mass inside the virial radius over 
the virial mass. Row 13: virial stellar fraction, i.e., stellar mass inside the virial radius over the virial mass. Row 14: stellar mass in solar masses. This is the 
stellar mass of the central galaxy. Row 15: effective stellar mass radius, i.e., half stellar mass radius in kiloparsecs. Row 16: stellar iron over hydrogen ratio. 
Mass weighted iron over hydrogen ratio for the dwarf stellar mass component. Row 17: total mass inside 500 parsec in solar masses. Row 18: dark matter 
mass inside 500 parsec in solar masses. Row 19: baryon mass inside 500 parsec in solai* masses. Row 20: gas mass inside 500 parsec in solar masses. Row 


^ simulation presented in 

Hopkins et al. 

2014 

as m09 

^ simulation presented in 

Hopkins et al. 

2014 

as mlO 


2 : 



Figure 4. Left panel: Star formation rates for the hydrodynamical runs presented in this work obtained using two different time bins: 10® yr (full red lines) 
and 10^ yr (dashed black lines). Right panel: Normalized cumulative SFR history for all the hydrodynamical runs presented in this work. Notice the difference 
of the star formation histories between the three high resolution runs. The vertical grey dashed line marks the reionization redshift assumed in the simulation. 
Simulated star formation histories are very similar to some observed ones (e.g.,|Skillman et al.|2014||Weisz et al.|2014^. See text for more details. 
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Figure 5. Evolution of the virial baryon fraction for all the hydrodynamical 
runs presented in this work. In all cases, gas is slowly expelled out of the 
halo by the stellar feedback. Notice how the Dwarf runs with the highest 
SFR at early times (when the halo is still growing) end up having a lower 
baryon fraction. The Ultrafaint is further affected by the UV background, 
which prevents accretion of gas after 2 : ~ 6. 


or virial velocity. The central potential depth increases only weakly 
over this time as the halo accretes material mostly on its outskirts. 
This low but constant SFR at low redshift is a key factor in shaping 
the final matter structure of the dwarf galaxy and will be discussed 
in detail in the next Section. 

Observations of the star formation rate of dwarf galaxies ^Tol 


stoy et al. 2009 Skillmanet al.|2014[|Weisz et al.|2014[|Brown et~ 

2014[ Co e et al.|2014 1 show a relatively high dispersion for a fixed 
stellar mass, but all the histories of our simulated galaxies seem re¬ 
alistic when compared with these data. In particular, our Ultrafaint 
galaxy is composed of uniformly old stars, as observed in real ultra¬ 
faint dwarf satellites of the Milky Wa5jBrown^et^J2014J(^erhaps 
making them fossils of reionization (|Ricotti & Gnedin|2005||Bovill| 
|& Ricotti|20rT) . 


Some insight into the extremely low efficiency of star forma¬ 
tion in all of these systems can be gained from examining the total 
baryon fraction vs. time within their associated virial radii. This is 
shown in Figure]^ Specifically, the virial baryon fraction (baryon 
mass divided by total mass inside the virial radius) begins declin¬ 
ing in the Ultrafaint run from the moment the UV background starts 
acting to reduce the amount of gas falling into the halo. This allows 
feedback to be more efficient in expelling the gas out of the halo po¬ 
tential. The Dwarf runs also begin to demonstrate a steady decline 
in their baryon fractions, but only after 2 ~ 2.5. This is the redshift 
when the halo itself stops growing (the end of the “rapid accretion 
phase”). Below that redshift, star formation feedback steadily acts 
to expel gas. The runs with the higher star formation rates have a 
lower baryon fraction, though all three of the Dwarf runs end up 
at a factor of ~ 10 — 20 below the cosmic mean (horizontal gray 
line). However, the stellar feedback in the DwarLearly simulation 
has managed to expel a larger fraction of material from halo, slow¬ 
ing down late time star formation. 

One intriguing result from Figure 5 is that the overall baryon 
fraction decreases steadily, without global jumps that are tightly 
linked to the star formation rate (which varies substantially over ~ 
100 Myr timescales). Instead, the baryons slowly “evaporate” out. 
This is in contrast with other studies ([Sawala et al.|20ir||Simpson| 


|et al.|20l'3) that show s harper jumps. However, the lower resolution 
used by Sawala et al. j20]^ could explain why their evolution is 
less smooth. ISimpson et al.| ( [20T^ used comparable resolution to 
this work but studied a lower mass system, 10®Mq, which could 
explain the much more drastic effect due to the UV background in 
the gas virial fraction that they found. 

In the next Section, we study how these differences in the star 
formation histories affect the matter distribution of the halo. 


4 DARK MATTER CONTENT AND STRUCTURE 

In Figure]^ we present the dark matter density profiles of the hy¬ 
drodynamical Ultrafaint (left) and Dwarf (right) rans compared 
with their equivalent collisionless run at 2 = 0. The grey bands 
mark the regions where the simulations are not fully converged ac¬ 
cording to the criterion of |Power et al.|j200T| computed for the dark 
matter only simulation. The Power radius is defined to be the radius 
where the two-body relaxation time, freiax, becomes shorter than 
the age of the universe to, where frelax is determined by the num¬ 
ber of particles and the average density of the enclosed region p. 
Specifically, Power et al. found that freiax < 0.6fo is the best crite¬ 
rion. [Elberrer^(|20^ have recently confirmed that this criterion 
is accurate using zoom simulations of collisionless dwarf halos at 
similar resolution to those we examine here. The vertical black dot¬ 
ted lines in Figure mark four times the dark matter gravitational 
softening used in the collisionless runs. We note that while radii 
larger than the Power radius should not suffer from two-body re¬ 
laxation, the smallest radius where results in hydrodynamical sim¬ 
ulations are converged may be (significantly) larger. 

The left panel of Figure|^shows results for the Ultrafaint sim¬ 
ulations. In this case, there is no sign of a decrease in the dark 
matter density in the hydrodynamical run; in fact the dark matter 
profile matches perfectly with the collisionless run. In the Dwarf 
runs (right panel), we observe that all hydrodynamical runs show 
varying levels of decrease of the inner dark matter density when 
compared with their equivalent collisionless run. In particular, the 
Dwarf date run has produced a fairly large (~ 1 kpc) constant den¬ 
sity core - this is exactly the behavior needed to help alleviate the 


Too Big to Fail Problem (see, e.g. |Elbert et al.||2014| |Governato| 
|et al.|[20r5) and that would be required to explain indications of 


cored profiles in low-mass galaxies in the Local Group 

Donato 

et al.|2009l|Salucci et al.|2012l|Walker & Penarrubia|2011 

lAmor- 

isco et al.|2014[|Burkert|2015|. 


One common way to quantify core formation in halos is to 
measure the log-slope of the density profile a at 1-2% of the 
virial radius. The Dwarf halo in the dark-matter only run has 
a = —1.58, while DwarLearly, DwarLmiddle, and DwarOate 
have a = —1.39, —0.88, and —0.27, respectively. The late- 
forming dwarf produces the shallowest profile and the largest core, 
while the early-forming dwarf produces the densest, cuspiest sys¬ 
tem. 

Over time, in all three dwarf runs, we have observed clear cor¬ 
relations between core formation and star-formation events. How¬ 
ever, at early times, as the halos continue to accrete matter and ex¬ 
perience central mergers, the cusps regrow regularly. During the 
early, rapid accretion phase, evolution in the density structure is 
fairly stochastic, with cores forming in response to blow-out events, 
and then becoming erased as cusps reform in response to mergers. 
Figurej^illustrates the formation of the dark matter core using two 
time steps. Shown are the dark matter density profiles of the Dwarf 
runs at 2 = 3.9 (left) and 2 = 2.2 (right). At 2 = 3.9, very lit- 


© 0000 RAS, MNRAS 000, 000-000 

















































Forged in FIRE 9 



Figure 6. Left: The dark matter density profile at 2 : = 0 for the collisionless (black line) and hydrodynamical (red line) runs of the 3 x 10^ Mq halo. The 
“collisionless” line has been converted to the effective dark matter density by accounting for the fact that a fraction of each particle is assumed to be 

baryonic in these nrns. The bottom panel shows the ratio between the two profiles. Right: The same for the Dwarf halo nins, where each hydrodynamical run 
is marked by a different style of red line. Grey shaded area marks the region below the convergence radius defined using |Power et al.]j2003t criteria for the 
colissionless run. The vertical black dotted line marks four times the dark matter gravitational softening used in the collisionless runs. Note that the Dwarf Jate 
run has produced a large (~ 1 kpc) constant-density core, while the Dwarf.early has a dark matter profile that is very similar to the dissipationless simulation 
for radii that are well converged. The dark matter in the hydrodynamic Ultrafaint run is identical to that of the dissipationless case. 



Figure 7. Time variation in density profiles. The dark matter density profile at z = 3.9 (left figure) and z = 2.2 (right figure) for the collisionless (black) 
and hydrodynamical (red lines) runs of the 1 X 10^® Mq halo. Bottom panel shows the ratio between the two profiles. The vertical black dotted line marks 
four times the dark matter gravitational softening used in the collisionless runs. Grey shaded area marks the region below the convergence radius defined using 
[Power et al.H2003) criteria for the colissionless run. Note that, at z = 2.2, Dwarfjate has a higher central density than Dwarf.early. Late-time star formation 
in Dwarfjate serves to reduce the dark matter halo’s density in the center by a factor of ~ 5 by z = 0 (see right panel of Fig. 6), while DwarLearly has little 
star formation subsequent to z = 2.2. Its density profile remains essentially unchanged from z = 2.2 to z = 0. 


tie star formation has occurred and the halo is experiencing very 
rapid growth and we see no decrease in core dark matter density 
compared to the collisionless run. However, at z = 2.2, there are 
some signs of a decrease in the central dark matter density in the 
hydro runs. Interestingly, Dwarfjate - which has the largest core 
at z = 0 - has the smallest core profile at z = 2.2. DwarLearly 
shows almost the opposite trend, owing to the fact that it has had 
more star formation by z = 2.2 than the later forming dwarf. 

To further explore the evolution of the dark matter density 
with time, Figure compares the cumulative star formation his¬ 
tory (dashed curves, normalized to unity at the present day) and the 
mass interior to radii of 0.3, 0.75, and 2 kpc relative to the colli¬ 
sionless run as a function of time for Dwarfjate (left panel) and 
DwarLearly (right panel). In both cases, the early phases of galaxy 


formation (z > 3) result in fluctuations in the inner mass profiles 
of these galaxies ( [Davis et al.|2014[ >. After z = 3, when the dark 
matter assembly of each halo is essentially complete, DwarLearly 
forms only a relatively small amount of stars. This results in at 
most a slight reduction in the inner dark matter mass (right panel). 
Dwarfjate, however, forms more than 50% of its stellar mass after 
z = 2. Most of the density reduction also occurs after this phase, 
pointing to a link between the final densities of these objects and 
their late-time star formation histories. This is consistent with ILa-l 
jporte & Penarrubia| ( |2015^ , who found that cusps can regrow after 
early core formation. 

Figurej^further illustrates the correlation between star forma¬ 
tion history and core formation, now with the early and late runs on 
the same plot, and using dimensional star formation histories rather 
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Figure 8. CoiTelation between star formation history and core formation for Dwarfjate (left) and Dwarf.early (right). The green dashed lines are the cumulative 
star formation histories for each run normalized to its value at 2 : = 0. The vertical line indicates the end-phase of the rapid-accretion epoch — feedback after 
this time is most capable of forming stable cores that are not regrown by subsequent mergers. The solid lines show the integrated mass within different radii 
(as labeled) divided by the same mass in the dissipationless runs - values less than unity indicate the formation of a lower density core than in the dark-matter 
only case. We see that in the Dwarfjate case (left), the core begins to form after the rapid accretion phase has ended (2 ~ 2) and keeps growing slowly down 
to redshift 2 = 0 as star formation continues. DwarLearly, however, forms most of its stars prior to the quiescent accretion phase and a large core never takes 
hold. 


than normalized ones. Specifically, the cumulative star formation 
histories of the DwarLearly (green dash) and Dwarfjate (red dash) 
runs are shown along with the evolution of the ratio of dark mat¬ 
ter enclosed within 0.3 kpc for the hydrodynamic compared to the 
dark-matter only runs (solid lines). It is clear that a higher star for¬ 
mation rate at late times (from 2 ~ 2) produces a bigger decrease 
in the central dark matter density. Notice that although this differ¬ 
ence in the star formation rate below z = 2 produces very different 
cores, the difference in the total amount of stars at 2 = 0 is minimal 
(see Figure|^and Table[^. In concordance with this picture, lower 
resolution runs, which have slightly higher star formation rates at 
low redshift than their high resolution counterparts, show bigger 
cores in their dark matter distribution (see Appendix [B). 

It is likely that the relationship between when the stars form 
and core formation is most important at this critical stellar mass / 
halo mass scale (M* ~ 2 — 3 x 10®Mq within ~ halos). 

Previous simulation efforts ( [Pontzen & Governato|2012||Di Cintio| 
|et al.|20l4) have found that dark matter cores are usually not cre¬ 
ated in galaxies with so few stars in halos below ~ IO^^Mq. We 
suggest that at this critical mass scale, where the energy from feed¬ 
back sources is just at the edge of that required for core formation, 
small variations in star formation histories can significantly alter 
the result. Indeed, in our general analysis of the dark matter prop¬ 
erties in all FIRE runs, we find a similar transition around 10^° Mq 
(C han et al, in preparation). 


4.1 Energy considerations 


Recently, there has been some discussion in the literature about the 
energy requirements for the formation of a core in a dwarf galaxy 
halo (see, for example Penarrubia et al.|[2012[ |Garrison-Kimmel| 


|et al.|20r3{ |Teyssier et al.|2013| and references therein) - specifi¬ 

cally, how many stars are required for there to be enough energy 
available to create a core? At first comparison, the fact that the 
Dwarfjate simulation was able to produce a sizable core with so 
few stars appears to be in contradiction to the results of Garrison- 


Kimmel et al. (2013), who suggested that cores this large are not en¬ 
ergetically possible. However, the host halo considered in Garrison- 
Kimmel et al. (2013) is more concentrated than the one we consider 
here. In order to explicitly check whether our results make sense en¬ 
ergetically we aim to compare the energy released in supemovae in 
our simulations to the difference in dark matter gravitational energy 
potential of the dwarf hydro runs and its collisionless version. 

The gravitational potential energy is defined as: 

y’^max 

U = —JttG / rM{r)p{r)dr (1) 

J ’’min 

We computed this value numerically directly from the simulation 
data. We considered /min = 0 (using rmin = Cdm or rmin = 
gives very similar results) and /max = 2 kpc. We used this max¬ 
imum radius because we are interested in the energy necessary to 
decrease the inner part of the density profile and from this point the 
dark matter profiles match almost exactly (see Figure]^. More im¬ 
portantly, at larger radii differences between the gravitational en¬ 
ergy potential can be significant just due to the exact position of 
the substructure between the collisionless and the hydrodynamic 
runs. Therefore this definition of potential energy for each run sets 
a lower energy limit on the amount of energy necessary to create a 
specific dark matter decrease in the inner part of a halo. We define 
AUdm as the difference between the potential energy of the hydro- 
dynamical run and the potential energy of the collisionless run. 

In order to obtain an estimate of the energy available from 
feedback we have considered the energy available from SNe using 
the parameters from our simulations: Etot ~ {Mata.i/m) f Ean, 
where Ean = 1 x 10®^ erg is the energy of one supernova, 
/ = 0.0037 is the fraction of stars more massive than 8 Mq 
for a ( |Kroupa|r2002^ IMF, and m = 0.4 M© is the mean stellar 
mass. The stellar mass of the central galaxy at 2 = 0 is between 
~ 2.3 — 2.8 X IO^Mq however from Figurej^we can see that the 
core starts to form below 2 ~ 2, therefore we have also consid¬ 
ered the stellar mass produced since this time until 2 = 0. We are 
considering just supernova energy but in principle, just taking into 
account the energy, the contribution from photoionization and radi- 
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Figure 9. The cumulative star formation history (dashed lines) and the dark 
matter mass ratio between the hydrodynamical run and the collisionless run 
at 0.3 kpc (full lines) for the early and late forming Dwarf runs. 


ation pressure could play a role in core creation. Although most of 
the energy in radiation just escapes the system, there is in principle 
~ 100 times more energy in radiation than in supernovas. Prelimi¬ 
nary tests done in this regard by changing the energy per supernova 
do not point towards a relevant role of these processes, at least as 
they are currently implemented in the code. However we leave a 
more careful analysis to future work. 

In Figure [T^ we plot the potential energy difference between 
the hydrodynamical runs and the collisionless run, A17dm, versus 
the total stellar mass (green symbols) and the stellar mass pro¬ 
duced since z = 2 down to 2 = 0 (red points) for all the different 
Dwarf runs. The supernova energy available from the stellar mass 
and the size of the cor^ linked with the difference in energy are 
also marked in the Figure. Grey lines stand for different efficien¬ 
cies of the supernova (1%, 5%, 10% and 30%), i,e. the energy from 
supernova affecting the dark matter. Our simulations indicate that 
the gravitational potential energy of the halo has changed by an 
amount consistent with ~ 5 — 10% of the SNe energy available 
for the Dwarf date run, while only ~ 1% of the total SNe energy 
has been effectively captured by the dark matter in the DwarLearly 
run, owing to the fact that many of these supemovae exploded dur¬ 
ing the period of rapid accretion, when cusps were reforming. We 
find similar results when considering the lower resolution runs, in¬ 
cluding efficiencies, as they have slightly bigger cores than their 
counterparts but a higher star formation rates at lower redshift (see 
Appendix[B). Overall, these results suggest that large cores induced 
by star formation feedback will never appear in galaxies smaller 
than ~ IO^Mq (in ~ 2 — 3 x IO^Mq or bigger halos), owing to 
the inefficiency of SNe energy coupling to the dark matter. 


5 SUMMARY AND CONCLUSIONS 

We have performed several high-resolution zoom-in hydrodynam¬ 
ical simulations of an ultrafaint galaxy halo (3 x 10® Mq) and a 
dwarf galaxy halo (1 x 10^® M©). Our simulations include all major 
sources of stellar feedback, implemented directly from stellar evo¬ 
lution calculations. Without parameter tuning, the code reproduces 


® We define the size of the core at the radius where the mass ratio between 
the hydrodynamical over the colissionless runs is 0.9 (see Figurej^. 


Figure 10. Energy considerations in the formation of the dark matter cores 
in a 10 ^®Mq halo. The x-axis show the total stellar mass for each of the 
dwarf runs (green symbols) and the amount of stellar mass formed between 
redshift 2 = 2 and 2 = 0 (red symbols). The upper axis show the energy 
expelled by SN from this stellar population. The y-axis show the difference 
in potential energy between the hydrodynamical run and the collisionless 
run. The right y-axis show the size of the core created due to this difference 
of energy. See text for details. 


a relation between galaxy stellar mass and halo mass that is consis¬ 
tent with observations. Moreover, we find that global properties of 
these simulated halos - including their characteristic sizes, metal- 
licities and gas contents - are well-matched to observed galaxies 
of similar stellar mass. These global properties describing the sim¬ 
ulated dwarfs are robust to changes in force and mass resolution. 
Furthermore, the feedback models and the outflows they generate 
are inherently multi-phase, matching observations. The predictive 
nature of our galaxy formation model is particularly important, as 
the model does not contain ad hoc numerical solutions adopted by 
other models, e.g. cooling shut-offs or prescribed wind properties, 
that contain adjustable parameters. The mass scale of our simulated 
dwarfs - M* ~ 2 x 10® Mq - is particularly relevant because pre¬ 
vious models able to generate cores have usually formed almost an 
order of magnitude more stars in such halos (Mvir = 10®° Mq). 
Such galaxies are too massive in terms of the number of stars given 
their halo mass and therefore cannot be typical, given observed 
galaxy counts around the Milky Way and generic predictions from 
ACDM simulations ( |Garrison-Kimmel et al.|[2014[ [Brook et al.| 
|20T4l l. 

Our models show a slow but continuous decrease of the bary- 
onic mass inside the virial radius after 2 ~ 6. The UV background, 
in concert with star formation feedback, plays a fundamental role 
in regulating star formation in low-mass systems and appears to be 
the driving factor in suppressing gas accretion in our ultrafaint run. 
However, for the halo masses studied in this work, the UV back¬ 
ground does not shut down star formation immediately because it 
is not efficient in heating the high density gas in the center of these 
halos. The simulated ultra-faint (Mvir = 3 x 10° M©) continues 
forming stars for ~ 2 Gyr following reionization (at which time it 
runs out of cold gas; such an object would be a counterpart in the 
field to known ultra-faint satellites of the Milky Way), while the 
more massive dwarfs continue to form stars to 2 = 0. This may 
indicate a transition from lower-mass objects that are incapable of 
acquiring cold gas after reionization to dwarfs at this mass scale 
that can continue to accrete fuel for subsequent star formation. 

We have also studied, in detail, the dark matter distribution of 
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these halos. The simulated dwarfs (M* ~ 2 x 10® M©) have a va¬ 
riety of density profiles, ranging from a small modification of the 
equivalent dark-matter-only simulation to a substantial (kpc-scale) 
core. The simulated ultra-faint galaxy (M* ~ 2x10'* M©) does not 
form enough stars to modify its dark matter halo at all, providing 
further support to the idea that there is a critical mass below which 
core formation caused by stellar feedback is energetically impossi- 
ble (see, e.g.,|Govemato et al.|2012]|Garrison-Kimmel et al.|2013[ 
|Madau et al.|2014| (. Our results indicate that stellar mass is not the 

only parameter in core creation, however. The creation of dark mat¬ 
ter cores is linked with late-time star formation properties, as only 
the system with significant late-time star formation forms a sizeable 
core. The galaxy that forms most of its stars at early time is able to 
create a core temporarily, but subsequent dark matter accretion and 
mergers and the lack of strong star formation erase this core, leav¬ 
ing a cuspier profile. The difference in density at 300 pc between 
these two extreme cases is a factor of « 4. A related point is that 
the formation of stable dark matter cores is a continuous process, 
not instantaneous, and that the creation of significant cores in dwarf 
galaxies does not appear to be an inevitable outcome in models with 
bursty star formation histories. 

A question that remains unclear is whether these cored sys¬ 
tems can avoid regenerating a density cusp once they merge with 
smaller, cuspier, haloes jLaporte & Pena rrubi a|2015l >. The late-time 
merger history of dwarfs can vary significantly (e.g., |Deason et al!] 
|2014[ >, meaning it is imperative to simulate a statistical sample of 
halos at a given mass to fully understand trends in core creation or 
cusp regrowth (Fitts et al., in preparation). It will also be imper¬ 
ative to test this scenario at different halo masses (Chan et al., in 
preparation), as many models predict a core formation efficiency 
that varies with the halo mass (e.g., |Di Cintio et al.|2014^ . We have 
not considered the effects stripping from ram pressure and tides, 
that may be important for some Milky Way subhalos ( [Read et ^ 
|2006|[Glotov et al.|2012[|Brooks & Zolotov|2014| l. However, the 
central prediction coming from our simulations is observationally 
testable: the presence of cores in galaxies with stellar masses of 
~ 10® — 10^ M© requires substantial late time star formation. 
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APPENDIX A: DARK MATTER PROPERTIES AND 
EVOLUTION IN THE COLISSIONLESS RUN 

To choose the specific dwarf galaxy halos to re-simulate we rely on 
colissionless simulations. We have taken into account two things: 
first we wanted them to be cheap in terms of cpu cost and second 
we wanted them to be representative of the dwarf galaxy halo pop¬ 
ulation. We point to |Onorbe et al.| ( [2014| l for a full description of 
the method. Here we just want to show how the properties of our 
selected halos compare with a realistic sample of dwarf galaxy ha¬ 
los. To generate this sample we mn a Lbox = 35 Mpc collisionless 
simulation (512^ simulations). Rigure Al show the spin (A), con¬ 
centration {VmaxIV-„ir), halo formation time (fso) and virial mass 
distributions for all the main halos in this simulation (so excluding 
subhalos) with virial masses between 3 x 10® Mq and 3 x 10^® 
Mq. The mass bin sample includes around ~ 15000 halos. Halo 
spin parameters were calculated using |Bullock et aL]p001| > defini¬ 
tion. In order to estimate the time of formation for each halo, we 
followed the approach described in |Wechsler et al.| ( |2002| ). We fit 
the halo accretion histories obtained from the merger trees to a ex¬ 
ponential form that depends on one parameter. The halo formation 
time tso is calculated at the time when the halo reached half of its 
total mass. The chosen parameters for our ultrafaint and dwarf ini¬ 
tial conditions are plotted as a white triangle and a white square 
respectively. The exact values can be found in Table[T] This Figure 
show that the re-simulated halos picked from our Li,ox = 7 Mpc 
box have very typical values of spin and concentration. The reason 
why our formation time is a bit lower than the standard value is a 
combination of three factors. First we preferred to avoid systems 
with late major mergers events which also helps to reduce the cpu 
cost of the simulation jOfiorbe et al.|2014) . The dwarf halo sample 
from a smaller box is biased towards smaller formation times so 
there were a smaller range of possible halos to pick which fulfill 
all our desired criteria. Finally, the circular velocity profile of the 
DwarLdm simulation, which can be found in [Elbert et al.| j2014[ > 
(left panel), shows that the halo has too high density to match the 
circular velocity observations of Local Field dwarf galaxies. This 
makes it a suitable candidate to study the Too Big To Fail problem. 

Figure fA^ shows the evolution of the dark matter mass profile 
for the Dwarf collisionless simulation. Each line shows the amount 
of dark matter mass contained inside a fixed physical radius. At 
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Figure Al. Selecting the sample. Color map shows the probability density 
function of the dwarf halo sample in the = 35 Mpc collisionless 

simulation (512® resolution). The white triangle and white square stand 
for the Ultrafaint.dm and DwarLdm runs respectively. Upper panel shows 
concentration versus halo spin. Lower panel shows the virial mass versus 
halo formation time. Specific values of these parameters can be found in 
Table0 


high redshift the halo shows a characteristic fast halo mass increase 
followed by a very shallow evolution at high redshift. Notice how 
the inner parts of the profile takes a bit more time to settle down. 
Below redshift 2 ; ~ 2.5 the inner part of the halo does not show 
any significant perturbation as there is no significant accretion or 
merger ( [Diemand et al.|2007l|Diemer et al.|20T^ . 


APPENDIX B: CONVERGENCE 

In this Section we present a convergence study that we have per¬ 
formed for the Dwarf galaxy halo. We have run a lower resolution 
version of all the Dwarf hydrodynamical runs discussed above. The 
only difference between the high and low resolution runs are the 
different particle masses and softenings used in the runs. The star 
formation density threshold, riaf, is also slightly different between 
the runs, 10 cm“® for low res and 100 cm“® for high res. All other 
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Figure A2. The evolution of the dark matter of the Dwarf collisionless run. 
Each line shows the dark matter mass contained inside a fixed physical ra¬ 
dius. 



code and physical parameters are exactly the same as for the high 
resolution runs. In Table |B1| all the relevant parameters of these 
lower resolution runs can be found. 

The different panels of Figure [Blj illustrate the differences be¬ 
tween the runs. The main difference that we found is that the low 
resolution runs have slightly higher stellar masses (upper left panel 
of Figure [bT) . This can be understood by looking at the SFR his¬ 
tories (upper right panel in Figure [BT| blue lines stand for low res¬ 
olution runs and red lines for the high resolution ones). The main 
difference observed between resolutions is the steeper slope of the 
cumulative star formation history at lower redshift. This produces 
higher stellar masses at 2 = 0 for the lower resolution runs. We 
think that this is because the minimum amount of star formation 
that is possible is set by the gas particle resolution. Therefore the 
minimum amount of star formation is higher in the lower resolution 
runs. It is remarkable that all galaxy trends with size and metallicity 
hold regarding of resolution, so the galaxies seems just move along 
these relations (lower left panel of Figure [BT) . We have also re¬ 
run our dwarf galaxy low resolution initial conditions using exactly 
the same code to check for pure stochastic differences. The scatter 
found in all the properties studied in this paper was similar to the 
one that we found when we change the feedback implementation 
and/or the softening values. These authors suspect that these differ¬ 
ences will decrease at higher resolution, though higher resolution 
runs will certainly be required in order to test this conjecture. 

Finally, concerning the core formation and energy consider¬ 
ations, low resolution runs also form a core which seems to be 
directly connected with its star formation rate at low redshifts 
(2 < 2). In general these cores are more prominent than their high 
resolution counterparts due to the higher star formation rate dis¬ 
cussed above. Remarkably when we plot the energy requirements 
to form these cores versus the amount of energy obtained from su¬ 
pernova feedback below 2 = 2 (lower right panel of Figure |BTJ, 
they lie in the same range of efficiencies as their high resolution 
counterparts. 

Although it is not possible to claim full convergence for our 
high resolution runs from these results, we think that they are at 
least quite encouraging and definitely an improvement from other 
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Figure Bl. Convergence tests. High (red) and low (blue) resolution simulations. Upper left: The stellar mass - halo mass relation. Upper right: Cumulative 
star formation history. Lower left: Metallicity versus stellar mass. Lower right: Energy considerations in the formation of the dark matter cores in a IO^^Mq 
halo. See text for details. 


approaches in which parameters of the suh-grid physics must he 
tuned at each resolution. 
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Parameter 

Dwarf_dm_lr 

Dwaif JateJr 

Dwarf_middle Jr 

DwarLearlyJr 


(Collisionless) 

(Hydro: Feed-M) 

(Hydro: Feed-M-soft) 

(Hydro: Feed-V) 

l)mdm 

1.21 X 10"^ 

1.01 X 10-* 

1.01 X lO-* 

1.01 X 10"* 

2) Cdm (pc) 

35 

35 

35 

35 

3) mbar 

- 

2.04 X 10^ 

2.04 X 10® 

2.04 X 10® 

4) c™ (pc) 

- 

2.0 

35 

2.0 

5) Mvir (Mq) 

9.48 X 10® 

7.60 X 10® 

7.60 X 10® 

7.46 X 10® 

6) Vmax (km/s) 

37.31 

32.79 

32.79 

33.56 

V Cvir (kpc) 

54.99 

51.08 

51.08 

50.77 

8) /bar 

- 

0.0166 

0.0171 

0.0137 

9) /gas 

- 

0.0160 

0.0165 

0.0121 

10 )/* 

- 

0.0006 

0.0006 

0.0016 

11)M* (Mq) 

- 

4.1 X 10® 

4.2 X 10® 

1.0 X lO’’ 

12) r*^2 (kpc) 

- 

0.783 

0.881 

1.311 

13) [Fe/H] 

- 

-1.493 

-1.468 

-1.450 

14) (Mq) 

3.212 X W 

1.285 X 10'^ 

1.014 X 10'^ 

4.088 X 10® 

15)Mf^*(M0) 

2.414 X 10^ 

7.366 X 10® 

7.132 X 10® 

3.703 X 10® 

16) Mb“(M0) 

- 

5.483 X 10® 

3.002 X 10® 

3.861 X 10® 

17) (Mq) 

- 

5.083 X 10® 

2.635 X 10® 

0.0 

(M©) 

- 

4.004 X 10® 

3.677 X 10® 

3.861 X 10® 


Table Bl. Simulations data for the low resolution convergence tests. First column stand for the different parameters studied for each simulation. In Columns 
2-9 results for the simulations presented in this work are shown. Row 1: dark matter particle mass in the high resolution region in solar masses. Row 2: fixed 
gravitational softening used for the dark matter paiticles in physical parsecs. Row 3: baryon particle mass in the high resolution region in solar masses. Row 
4: minimum baryonic force softening in parsecs (minimum SPH smoothing lengths are comparable or smaller). Recall, force softenings are adaptive (mass 
resolution is fixed). Row 5: virial mass in solar masses defined at the overdensity at which the spherical top hat model predicts virialization jBryan & Norman] 
|1998) . Row 6: maximum circular velocity in km/s. Row 7: virial radius in kiloparsecs. Row 8: virial baryon fraction, i.e., baryon mass inside the virial radius 
over the virial mass. Row 9: virial gas fraction, i.e., gas mass inside the virial radius over the virial mass. Row 10: virial stellar fraction, i.e., stellar mass 
inside the virial radius over the virial mass. Row 11: stellar mass in solar masses. This is the stellar mass of the central galaxy. Row 12: effective stellar mass 
radius, i.e., half stellar mass radius in kiloparsecs. Row 13: stellar iron over hydrogen ratio. Mass weighted iron over hydrogen ratio for the dwarf stellar mass 
component. Row 14: total mass inside 500 parsec in solar masses. Row 15: dai'k matter mass inside 500 parsec in solar masses. Row 16: baryon mass inside 
500 parsec in solar masses. Row 17: gas mass inside 500 parsec in solar masses. Row 18: stellar mass inside 500 parsec in solar masses. 
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